STOCHASTIC EULERIAN-LAGRANGIAN METHODS FOR 

FLUID-STRUCTURE INTERACTIONS WITH THERMAL 
FLUCTUATIONS AND SHEAR BOUNDARY CONDITIONS 

PAUL J. ATZBERGER * 

Abstract. A computational approach is introduced for the study of the rheological properties of 
complex fluids and soft materials. The approach allows for a consistent treatment of microstructure 
elastic mechanics, hydrodynamic coupling, thermal fluctuations, and externally driven shear flows. 
A mixed description in terms of Eulerian and Lagrangian reference frames is used for the physical 
system. Microstructure configurations are represented in a Lagrangian reference frame. Conserved 
quantities, such as momentum of the fluid and microstructures, are represented in an Eulerian 
reference frame. The mathematical formalism couples these different descriptions using general 
operators subject to consistency conditions. Thermal fluctuations are taken into account in the 
formalism by stochastic driving fields introduced in accordance with the principles of statistical 
mechanics. To study the rheological responses of materials subject to shear, generalized periodic 
boundary conditions are developed where periodic images are shifted relative to the unit cell to 
induce shear. Stochastic numerical methods are developed for the formalism. As a demonstration 
of the methods, results are presented for the shear responses of a polymeric fluid, lipid vesicle fluid, 
and a gel-like material. 
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1. Introduction. Soft materials and complex fluids are comprised of microstruc- 
tures which have mechanics and interactions characterized by energy scales com- 
parable to thermal energy. This feature results in interesting bulk material prop- 
erties and phenomena which often depend sensitively on temperature and applied 
stresses 0; d; [l2| ■ Example materials include colloidal suspensions, foams, polymeric 
fluids, surfactant solutions, lipid vesicles, and gels [nl E; [H; [H; H; 32; 33]. Mi- 



crostructures of such materials include flexible filaments, bubbles, colloidal particles, 
lipid chains, and polymers. These microstructures are typically surrounded by a 
solvating fluid which further mediates interactions through solvation shells [26c l29l ] 
and hydrodynamic coupling [B|; 0; EH- I n addition, given the energy scales of the 
microstructure mechanics and interactions, thermal fluctuations often play an impor- 
tant role both in microstructure organization and kinetics [7|; 1 12c [32| . A fundamental 
challenge in the study of soft materials is to relate bulk material properties to mi- 
crostructure mechanics, interactions, and kinetics. 

For the study of soft materials we introduce a modeling and simulation approach 
which consistently accounts for microstructure elastic mechanics, hydrodynamic cou- 
pling, and thermal fluctuations. The modeling approach is based on a mixed Eulerian 
and Lagrangian description. The microstructure configurations are modeled in a La- 
grangian reference frame, while an Eulerian reference frame is used to account for 
conserved quantities, such as momentum, of the system. When coupling these dis- 
parate descriptions an important issue is to formulate methods which do not introduce 
artifacts into the conservation laws, such as artificial dissipation of energy or loss of 
momentum. These properties are especially important when introducing stochastic 
driving fields to account for thermal fluctuations. We discuss a general approach for 
developing such coupling schemes, focusing primarily on one such realization referred 
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to as the Stochastic Immersed Boundary Method [H; [3l[ . 

To facilitate studies of the rheological properties of soft materials we introduce 
methods to account for externally driven shear flows. To account for shearing of the 
material, we generalize the usual periodic boundary conditions so that periodic images 
are shifted relative to the unit cell to induce shear. Our approach is based on bound- 
ary conditions introduced for Molecular Dynamics simulations which are referred to 
as Lees-Edwards boundary conditions [27| . These conditions present a number of 
challenges in the context of numerically solving the hydrodynamic equations. We 
develop numerical methods which utilize jump conditions in the velocity field at do- 
main boundaries and utilize a change of variable to facilitate handling of the shifted 
boundaries. Further issues arise when accounting for the thermal fluctuations. For the 
introduced discretizations we develop stochastic driving fields which yield stochastic 
numerical methods which are consistent with the principles of statistical mechanics. 

We consider primarily two physical regimes. In the first, the relaxation dynamics 
of the fluid modes is explicitly resolved. In the second, the fluid modes are treated 
as as having relaxed to a quasi-steady-state distribution. In the first regime we de- 
velop efficient stochastic numerical methods for the generation of the corresponding 
fluctuating fields of the fluid. In the second, we develop efficient stochastic numerical 
methods which account for the correlated stochastic driving fields which account for 
the effective thermal fluctuations which drive the microstructure dynamics. 

As a demonstration of the proposed stochastic numerical methods, simulations 
are performed for specific systems. These include studying the shear responses of (i) 
a polymeric fluid, (ii) a vesicle fluid, and (iii) a gel-like material. To relate microstruc- 
ture interactions and kinetics to bulk material properties we develop estimators for an 
associated macroscopic stress tensor. The estimators take into account the n-body in- 
teractions in the microstructure mechanics and the generalized boundary conditions. 
For the polymeric fluid, this notion of stress is used to investigate the dependence of 
the shear viscosity and normal stresses on the rate of shear. The vesicle fluid is sub- 
ject to oscillating shear and simulations are preformed to characterize the frequency 
response in terms of the elastic storage modulus and viscous loss modulus over a wide 
range of frequencies. As a further demonstration of the methods, the time dependent 
shear viscosity of a gel- like material is studied through simulations. 

The ability to simulate explicitly the microstructure mechanics, hydrodynamic 
coupling, and thermal fluctuations provides an important link between bulk mate- 
rial properties and phenomena at the level of the microstructures. The presented 
framework and related stochastic numerical methods are expected to be applicable in 
the modeling and simulation of a wide variety of soft materials and complex fluids. 
The general approaches introduced for coupling the Eulerian and Lagrangian descrip- 
tions and for the incorporation of thermal fluctuations are expected to allow for the 
development of many different types of Stochastic Eulerian Lagrangian Methods. 

2. Stochastic Eulerian-Lagrangian Modeling Approach. We use a mixed 
Eulerian-Lagrangian description. The conserved quantities of the entire system in- 
cluding both the fluid and microstructures will be accounted for in an Eulerian refer- 
ence frame. The microstructure configurations will be accounted for in a Lagrangian 
reference frame, see Figure 12.11 To introduce the basic approach and simplify the 
presentation we consider here only a rather special case. We shall assume the solvent 
fluid is an incompressible Newtonian fluid of constant density and the microstructures 
are density matched with the fluid. In this case, the primary conserved quantity of 
interest is the local momentum of the material. The basic framework and principles 
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Fig. 2.1. ^4 description of the physical system is used which combines Eulerian and Lagrangian 
reference frames. The configuration of the micro structures are described using a Lagrangian refer- 
ence frame, shown on the left. Microstructures represented in the Lagrangian frame may be solid 
bodies, membrane structures, polymeric structures, or point particles. The conserved guantities, 
such as the local momentum, mass, or energy, are described in an Eulerian reference frame, shown 
on the right. The mapping X(q) relates the Lagrangian reference frame to the Eulerian reference 
frame. 



that will be presented are more generally applicable allowing for additional conserved 
quantities to be taken into account, such as the local mass density and energy [IBj . A 
more abstract and general presentation of the formalism will be the focus of another 
paper. 

The basic Eulerian-Lagrangian formalism describing the state of the fluid and 
microstructures is given by the following equations 

(2.1) DP m f) = v ' °" (x ' t] + A(x ' l) + A(x ' t] + g(x ' l) 

(2.2) ^l^ = r(q,t)+ 7 (q,t)+Z(q,t). 

The p accounts for the momentum of the material occupying location x and Dp/ Dt 
denotes the material derivative. The X(q, t) denotes the configuration of the material 
at time t and parameterized by q. The local material stress is denoted by cr = er[p, X]. 
We use the convention that cr accounts only for the dissipative stress contributions 
in the system. The operators A, T couple the Eulerian and Lagrangian descriptions 
of the state of the material. The operator A = A[X] accounts for momentum gained 
or lost locally in the system as the material deforms from non-dissipative stresses and 
body forces. The operator T = T[p, X] determines the rate of deformation of the 
material from the momentum of the system. The A = A[X,p] and 7 = 7[X, p] are 
Lagrange multipliers associated with time-independent kinematic constraints imposed 
on the system, such as rigidity of a body or incompressibility. Thermal fluctuations 
are taken into account through the stochastic fields g and Z. 
We consider systems where the total energy is given by 

(2.3) J B[p,X] = |^|u(x)| 2 dx + $[X], 

where u(x) = J o " 1 p(x) is the velocity of the material at location x, po is the constant 
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mass density of the material, and <£> is the potential energy for a given configuration. 
The force associated with this energy is denoted by F = — £$/(5X. 

For the operators which couple the Eulerian and Lagrangian descriptions to be 
physically consistent, the following should hold: (i) the coupling operators should not 
introduce any loss or gain of energy, (ii) momentum should only change through forces 
acting within the system. More precisely, these conditions require 

(2.4) J F(q) • F(q)dq = J [p„ x p(x)] ■ A[F](x)dx 

(2.5) ^A[F](x)dx = J P(q)dq. 

The conditions are required to hold for any realization of X, p, and F. The condi- 
tion 12.41 ensures the coupling operators conserve energy. The condition 12.51 ensures 
that in the absence of constraints the total momentum change of the system is equal 
to the total force acting on the system. 

In the notation, we find it convenient to write the operator A as explicitly de- 
pending on both X and F, which for conservative forces is technically redundant. To 
simplify the discussion, it has been assumed that the stress contributions denoted by 
<r are entirely dissipative and that there is no net in-flux of momentum from boundary 
stresses J dQ <r(x) • ndx = 0. 

With these conditions satisfied by the coupling operators, we discuss how to 
account for thermal fluctuations using the stochastic fields g and Z. It is convenient 
when accounting for thermal fluctuations to introduce coupling operators so that all 
configurations X are equally probable at statistical steady-state, when the $[X] = 0. 
It can be shown that this corresponds to dynamics determined by the constraints 
and coupling operators which introduces an incompressible flow on phase space. The 
requirement of an incompressibile flow on phase space can be expressed as 

(2.6) J — (x,x)dx + J — (x,x)dx + J _(q,q)dq + J ^L(q,q)dq = 0. 

This condition can be relaxed to allow for more general choices of coordinates, coupling 
operators, and constraints. If this condition is not satisfied a more general treatment 
of the thermal fluctuations is required to take into account in the invariant distribution 
the local compression or dilation of volume under the phase space flow [39| . 

To simplify the discussion, we assume that the dissipative processes can be ac- 
counted for by a negative definite self-adjoint linear operator C in p, so that Ver = Cp, 
and that conditions 1 2 . 41 - 12"1)1 are satisfied. With these assumptions, the thermal fluc- 
tuations can be accounted for using for g and Z Gaussian stochastic fields which 
are mean zero and 5-correlated in time [T§; |30| . The main issue then becomes to 
determine an appropriate spatial covariance structure for these stochastic fields. By 
requiring that the Boltzmann distribution be invariant under the stochastic dynamics 
of equations 12.11 - [2~2l it is required that Z = 0, and that 

(2.7) G(x, t, y, a) = ((g(x, t))(g(y, s)f) = -2k B T Po S(t - s)£5(x - y), 
see Appendix 1X1 

Similar formulations as equations 12.11 - 12.21 with g = 0, Z = 0, are the start- 
ing point for the derivation of a wide variety of computational approaches used for 
systems in which fluids interact with rigid or elastic bodies. These include Arbi- 
trary Lagrangian- Eulerian Methods (ALE) [Hill, Fluctuat ing Immersed Material 
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(FIMAT) Dynamics [301, Immersed Finite Element Methods (IFEM) [H |43|, and 
Immersed Boundary Methods (IBM) 0; [jnj]. The approaches we introduce allow 
for the incorporation and simulation of thermal fluctuations in such methods, which 
collectively we refer to as Stochastic Eulerian-Lagrangian Methods (SELMs). 

3. Semi-Discretization of the Momentum Equations, Microstructures, 
and the Eulerian-Lagrangian Coupling. We now consider semi-discretizations 
of the SELM equations. The momentum equations will be spatially discretized on 
a uniform mesh. The p m will denote the momentum at the mesh site indexed by 
m = (mi, 77i 2 , 7713) and the composite vector of such values will be denoted by p. The 
deformation state which describes the microstructure configurations will be discretized 
using a finite number of degrees of freedom denoted by indexed by j = 0, 1, . . . , M 
and the composite vector denoted by X. As an energy for this discretized system we 
use 

(3.1) £[p,X] =^i Po - 1 | Pm | 2 A a; d + $(X) 

m 

where Ax is the mesh spacing and d is the number of dimensions. The semi-discretization 
in space of the momentum and configuration equations can be expressed as 

(3.2) ^E = Lp + A + A + g 

Dt 

(3.3, 2*£W =rM+7M 

where D / Dt and L denote respectively the spatially discretized approximation of the 
material derivative and C. The p, A, A, g denote the composite vector of values on 
the mesh and X^, r^, 7^ denote values associated with the j th configurational 
degree of freedom. We assume the discrete dissipative operator is symmetric L = L T 
and negative semi-definite. For the coupling operators of the discretized equations 
the corresponding consistency conditions 12.41 - l2~5l are 

(3.4) T[p] T F = Pq 1 p T A[F]Aa; d 

(3.5) J2 A ^™ Axd = Y. FU] - 

m j 

The superscript T denotes the matrix transpose. The first condition ensures for the 
discretized system that the coupling preserves the energy and the second that changes 
in momentum only occur from forces acting within the system. The phase space 
incompressibility condition corresponding to equation 12.61 in the discretized setting 
becomes 

(3.6) V p -(A + A)+Vx-(r + 7 ) = 0. 

This condition ensures the uniform distribution for the configurations X is invari- 
ant under the stochastic dynamics of equation 13.2113.31 when the potential energy is 
constant, i.e. $ = 0. When A and T are linear operators the energy condition 12.41 
amounts to the coupling operators being adjoints (up to a scalar), 

(3.7) T = A T Po Ax d . 
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Provided conditions 13.41 - I3T6I are satisfied, the thermal fluctuations can be taken 
into account using a Gaussian stochastic field on the lattice, without any direct 
stochastic forcing required in the microstructure equations 12.21 If these conditions 
are violated by the discretization the numerical approximation may introduce artifi- 
cial loss or gain of energy or momentum in the system. In order to be consistent with 
the fluctuation-dissipation principle of statistical mechanics such discretizations would 
require additional sources of stochastic forcing to obtain the appropriate Boltzmann 
ensemble. 

An appropriate covariance structure for the stochastic driving field for discretiza- 
tions satisfying conditions 13.41 - 13.61 can be determined by requiring invariance of the 
Boltzmann distribution under the stochastic dynamics of equation 13.21 - 13.31 This 
yields for the semi-discrete system, see Appendix [A] 

(3.8) G = (gg T ) = -2LC. 

The covariance of the equilibrium fluctuations is given by the entries 

(3.9) C = P -^I, 

v ; Ax d 

where I is the identity matrix, see Appendix IA1 

One such realization of this SELM approach is the Stochastic Immersed Boundary 
Method developed in 0]. In this case the coupling operators are given by 

M 

(3.10) A /B F = £)F0(X(t))$ a (x m - X0(t)) 

3=1 

(3.11) [Ti B n] m = ^ a (x m - X^(t))u m (t)Ax d . 

m 

where u = Pq 1 P, and S a is a special kernel function approximating the Dirac 5- 
function, see Appendix [Cl The dynamics are subject to the constraint that the fluid is 
incompressible V-u = 0. For the semi-discretized SIB method of [H the conditions [3T] 
and 13.51 can be readily verified to hold exactly. However, the condition 13.61 only 
approximately holds and is exact only in the continuum limit. As Ax — > we have 

(3.12) V xW ■ r = J2 - V(5 a(x m - XM)umA/ - f -WS a (y - X^)u(y)dy 

= J £ a (y-XM)V-u(y)dy = 0. 

In the last line we used that the fluid is incompressible V • u = 0. The divergence of 
the terms A, A, 7, are zero in this case. In practice, the method still yields reasonable 
results since the exhibited fluctuations deviate from Boltzmann statistics only up to 
the discretization error 0; Q ■ 

4. Soft Materials and Complex Fluids Subject to Shear. We now discuss 
how the SELM approach can be used for the study of rheological properties of soft 
materials and complex fluids. We then discuss specific stochastic numerical methods 
for performing simulations in practice. For the type of materials we consider, it will be 
assumed that the solvent hydrodynamics is described well by the constitutive laws of 
Newtonian fluids in the physical regime where the Reynolds number is small. We also 
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assume that the explicitly represented microstructures occupy only a relatively small 
volume fraction and are effectively density matched with the solvent fluid. In this 
regime, the momentum of the system will be accounted for using the time-dependent 
Stokes equations 

(4.1) p -^ =/iAu- Vp + A + g 

(4.2) V • u = 

where u(x, t) = p^" 1 p(x, t) is the local velocity of the fluid body at x in the Eulerian 
reference frame, po is the fluid density, p is the dynamic viscosity, and p is the pres- 
sure. This corresponds to the dissipative stress a — [i (Vu + Vu T ) and A = — X7p 
in equation 12.11 While the Reynolds number is small, the partial time derivative is 
retained in the Stokes flow in equation 14.11 since the thermal fluctuations introduce 
small characteristic time scales into the dynamics. 

To introduce shear we generalize the usual periodic boundary conditions. Our 
basic approach is motivated by the molecular dynamics methods introduced by Lees- 
Edwards [3; 17; 27 1 . In this work, molecules in the base unit cell have modified 



interactions with molecules in periodic images. To simulate a bulk material undergo- 
ing a shear deformation at a given rate, the periodic images are treated as shifting in 
time relative to the unit cell, see Figure 14.11 This has the effect of modifying both 
the location of periodic images of molecules and their assigned velocities. This has 
some advantages over other approaches, where an affine-like deformation is imposed 



on the entire material body [18|; [24j; [40| • In contrast, for the Lees-Edwards approach 
the shear deformation is only imposed at the boundaries allowing within the unit cell 
for the molecular interactions to determine the form of the shear response. 

Motivated by this molecular dynamics condition we develop a corresponding 
methodology for the SELM approach. For momentum accounted for by the time- 
dependent Stokes equations we introduce the following generalized periodic boundary 
conditions 

(4.3) u(x, y, L, t) = u(x - vt, y, 0, t) + ve x . 

For concreteness we consider the case where a shear is imposed in the z-direction 
giving rise to velocities in the x-direction. The L is the side length of the periodic cell 
in the z-direction, v = Lj is the velocity of the top face of the unit cell relative to 
the bottom face, 7 denotes the rate of shear deformation, and ej is the standard unit 
vector in the j th direction. The interactions between microstructures of the system 
can be readily handled in the same manner as in the molecular dynamics simulation. 
This is done by shifting the location of any microstructure of a periodic image involved 
in an interaction. 

While conceptually straight-forward, these boundary conditions present signifi- 
cant challenges in practice for the numerical discretization of the momentum equa- 
tions. The conditions introduce both a jump discontinuity at periodic boundaries 
and a shift which potentially leads to misalignment of discretization nodes at the do- 
main boundaries, see Figure [4~T1 For commonly employed approaches such as spectral 
Fourier methods the jump discontinuity results in a degradation of accuracy through 
the resulting Gibbs' phenomena [21] . For uniform finite difference methods on the 
unit cell the mesh misalignment requires modified stencils or interpolations at the do- 
main boundary. When incorporating stochastic driving fields to account for thermal 
fluctuations these issues are further compounded. 
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Fig. 4.1. Discretization Mesh with Lees- Edwards Boundary Conditions Initial Formulation. 
The boundary conditions induce shear by shifting the periodic images of the unit cell, shown on the 
left and middle. For uniform discretizations of the unit cell this presents challenges since the mesh 
becomes misaligned at boundaries between the unit cell and the periodic images, shown on the right. 



To address these issues, we develop discretization methods which utilize a moving 
coordinate frame which deforms with the unit cell, see Figure 14.21 Let the velocity 
field in this frame be denoted by w(q, t) := u(0(q, i), t), where q = (q\, q 2 , q^) param- 
eterizes the deformed unit cell. Let </>(q, t) = (qi + q 3 A /t,q 2 ,q 3 ) denote the map from 
the moving coordinate frame to the fixed Eulcrian coordinate frame x = 0(q). The 
time-dependent Stokes equations in the deforming coordinate frame become 

(4.4) — — = Po V [e d - S d , 3 jte x f V 2 w^ [e d - S d , 3 jte x ] — Vp + F + J 

(4.5) V • w - ej Vw e x A /t = K 

where q — (qi, q 2 , q$) parameterizes the deformed unit cell, 7 denotes the rate of the 
shear deformation, the standard basis vector in the i direction with i G {x,y,z}. 
In the notation the parenthesized superscript denotes a vector component and 5k,e 
denotes the Kronecker (^-function. We also use the notational convention 



(4.6) 

dqidq-j 

(4-7) ^-^q-- 

In the equations, the terms J,K are introduced to account for the jump introduced 
by the boundary conditions 14.31 This allows in the new coordinate frame for use of 
the usual periodic boundary conditions 

(4.8) 

We now discuss a discretization for equations 14.41 and 14.51 and the corresponding 
source terms J, K. The following central finite difference approximations will be used 

,. n , dw( d ) w( d )(q + e l )-w( d )(q-e l ) 

(4.y) 



w(q 1 ,q 2 ,L,t) = yv(q 1 ,q 2 ,0,t). 
retization for equations 14.41 and 



dqt 2Ax 
02 w (d) w ( d )(q + e i +ej) - w( d )(q- 



^ 4 ' 10 ^ dqidq, ' 4Ax 2 
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Fig. 4.2. Discretization Mesh with Lees-Edwards Boundary Conditions Using a Moving Co- 
ordinate Frame. By discretizing the momentum equations in a moving coordinate frame a uniform 
discretization is obtained in which the mesh of the unit cell aligns with the mesh of the periodic 
images, shown on the left, middle, and right. The definition of the unit cell is changed from a cube 
to a sheared parallelepiped, shown on the far right on the bottom. 



(4.11) 



w^(q + ej-ej)- w% - e. t - ej) . 

4Ax 2 ' 1 * 3 

g2 w (d) w W( q + ei ) - 2wW(q) + w^fq - e 4 ) 
dqf Ax 2 



These approximations are substituted into equations 14.6114.71 to approximate the op- 
erators in equation 14.41 14.51 

We remark that the moving coordinate frame makes the description of the mo- 
mentum field have some features of a Lagrangian frame of reference. We none-the-less 
retain the Eulerian terminology treating this distinction loosely since the deformation 
corresponds to a somewhat arbitrary coordinate frame introduced for numerical conve- 
nience and does not directly follow from the details of the fluid flow. Our discretization 
approach shares features with Arbitrary Eulerian Lagrangian (ALE) Methods 0G5|. 



An important issue when using such deforming reference frames is that the dis 



cretization stencils may become excessively distorted [13J; iLM ■ We avoid this issue by 
exploiting the periodic symmetry of the system in the x and y directions. Let the 
displacement in the x-direction of the top of the unit cell relative to the bottom of the 
unit cell be denoted by s. For shear rate 7 and cell size L the displacement at time t is 
given by s = Ljt. The periodicity in the x and y directions has the consequence that 
for any coordinate frame with s > L there is another coordinate frame with s < L 
which has aligned mesh sites, see Figure 14.21 By adopting the convention that the 
coordinate frame with s < L is always used when evaluating stencils the distortion is 
controlled. 

To obtain approximations for the source terms J, K the discretization stencils 
are applied at the shear boundaries of the unit cell. For any stencil weights involv- 
ing values at mesh sites which cross the boundary the modified image value is used 
w m ± jL. The contributions of the stencil weights multiplied by ±jL are collected 
over all boundary mesh sites to obtain the source terms J, K. This allows for the 
usual finite difference stencils to be used on the unit cell with regular periodic bound- 
ary conditions. When including the source terms this gives the equivalent result of 
imposing the jump boundary condition 14.31 This formulation has a number of ad- 
vantages when numerically solving the discretized equations and when introducing 
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thermal fluctuations. 

The Stokes equations 14.41 - 14.51 discretized in this manner on a uniform periodic 
mesh can be expressed as 

(4.12) ^ =L(t)w + F + J + g 

(4.13) D(t)w = K 

where L(t) denotes the finite difference operator approximating the Laplacian in the 
moving coordinate frame, and D{t) the approximation of the Divergence operator 
in the moving frame, see equation 14.51 The discretization L{t) can be shown to be 
symmetric and negative semi-definite for each t. 

An important property of L(i),D(t) is that for each time t the corresponding 
stencils are translation invariant with respect to lattice shifts of the mesh. This has 
the important consequence that the matrix representations are circulant and therefore 



diagonalizable by Fast Fourier Transforms 37]. As a result, the incompressibility 
constraint can be handled using FFTs to obtain an exact projection method 0. This 
allows for the discretized approximation of the Stokes equations to be expressed as 

(4.14) = p(t) [L(t)w + F + J] + g 

where p(t) is the operator which projects to the null space of D(t). The incompress- 
ibility condition is then satisfied for all time provided D(0) ■ w(0) = K. 

We discuss stochastic numerical methods for two particular physical regimes: (i) 
the relaxation of the hydrodynamic modes of the system is resolved explicitly, (ii) for 
the current configuration X the hydrodynamic modes are treated as having relaxed 
to statistical steady-state. We remark that the case of resolving the hydrodynamic 
relaxation of the system is amenable to stochastic numerical methods similar to those 
introduced in We discuss this case only briefly and focus primarily on the newly 
introduced stochastic numerical methods for handling the second case. 

5. Regime I : Resolution of Hydrodynamic Relaxation. We now discuss 
in practice how the stochastic fields may be generated in the regime where the relax- 
ation of the hydrodynamic modes is resolved explicitly. For this purpose we express 
equation 14. 141 in differential form 

(5.1) dw = p(t) [L(t)w + F + J]dt + QdB t . 

The QdH t denotes the stochastic driving field accounting for thermal fluctuations cor- 
responding to g and Bf € E 3JV denotes the composite vector of a standard Brownian 
motion process at each of the mesh sites. Throughout our discussion the stochastic 
differential equations will be given the Ito interpretation (l9| . 

Using (Qd& t dR>t Q T ) = QQ T dt = Gdt, we see that Q denotes a matrix square- 
root of the covariance of the stochastic driving field G = QQ T . Given the discretiza- 
tions introduced in Section |4j the dissipative operator L(t) depends on time, which 
requires, see Appendix IB") 

(5.2) G(t) = -2p{t)L{t)C. 



This has the consequence that the covariance of the stochastic driving field is time 
dependent. 
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In the discretized system the numerical stencils dependent on time, However, since 
the shear deformation is volume preserving the discretized summation introduced to 
model the kinetic energy of the discrete system in equation 13.11 evaluated in the 
deformed coordinates is in fact not dependent on time. Since we made this choice, 
we have the important consequence that the Boltzmann equilibrium fluctuations of 
the velocity field w associated with this energy are stationary (independent of time) . 
In other words, the covariance C of the equilibrium fluctuations on the discretized 
lattice for the energy given in equation 13. II is independent of time, C{t) = C(0). This 
holds even though the underlying discretization and corresponding operators L(t), 
p(t), G(t) depend on time. 

To obtain an explicit form for G(t) we need to compute C taking into account the 
incompressibility constraint 14.131 The equilibrium covariance under these constraints 
is given by 

(5-3) C=*Wl. 
y 1 3 p Ax d 

The factor 2/3 arises from application in Fourier space of the projection operator 
which equivalently enforces the incompressibility. The factor po appears in the de- 
nominator since the velocity w is considered, instead of the momentum p = po~w. 
The notation for the stochastic driving field g is used loosely when switching between 
the momentum and velocity equations. 

The time dependent covariance structure of the stochastic driving field g in equa- 
tion [5TTJ is of the form G(t) = — 2p(t)L(t)C. An important issue is whether this will 
indeed yield a consistent treatment of the thermal fluctuations so that the resulting 
stochastic dynamical system has the required equilibrium fluctuations. We establish 
a Fluctuation-Dissipation principle for such time dependent systems in Appendix [Bj 

6. Generating the Stochastic Driving Field I. In order for equation 15.21 
and 15.31 to be useful in practice, we must have efficient methods by which to generate 
the stochastic driving fields with the required covariance structure. A significant 
challenge in practice is to generate efficiently the Gaussian stochastic driving field 
with the required covariance structure G(t). A commonly used approach is to generate 
a variate with uncorrelated standard Gaussian components £ and set g = Q(t)$, for 
an appropriately chosen matrix Q{t). The resulting variate g then has covariance 
(gg T ) = Q(t)(^ T )Q{t) T = Q{t)Q(t) T = G(t), with a proper choice of Q(t). 

However, to carry this out in practice encounters two challenges: (i) given G(t) 
the factor Q(t) must be determined, (ii) the matrix- vector multiplication Q(t)£ must 
be carried out. For (i) the Cholesky algorithm is typically used with a computational 
cost of 0(N 3 ), where N is the number of components of g. For (ii) the resulting 
factors Q(t) are generally not sparse, which when generating each variate incurs a 
computational cost of 0(N 2 ). To get a sense of the costs, for a three dimensional 
mesh, the number of components of g scales cubically as N = (£/Ax) 3 , where I is 
the domain size and Arc is the mesh resolution. The associated costs for generating 
the variates using this approach even for moderate spatial resolutions is prohibitively 
expensive. 

To obtain a more efficient computational method we use specific features of the 
discretization introduced in Section [4j One useful feature of the discretization we 
use is that the equilibrium covariance matrix is proportional to the identity matrix 
C = al with a — 2ksT /3poAx d . This allows equation 15.21 to be expressed as 



(6.1) 



G(t) = -2ap(t)L(t). 
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We also use the following specific properties of the operators C, L(t), and p(t) obtained 
from the discretization. The first is that each of the operators corresponds to use of 
numerical stencils which are translation invariant on the mesh. This has the important 
consequence that all of these operators are diagonalizable in the Fourier basis. This 
has the further important consequence that all of these operators commute. The 
second is that p is an exact projection operator, so that p 2 = p and p = p T . Finally, 
we use that the discrete approximation of the Laplacian is symmetric negative semi- 
definite so that it can be factored as L(t) — —U(t)U T (t) for some matrix U(t). 

By using these properties of the operators we can express the covariance of the 
stochastic driving field as 



(6.2) G(t) = (V2^pU(t)) (\/2^pt/(*)) 



r 



In this form the required matrix square-root is readily obtained as Q(t) — \f2apU (t) . 
We remark this is different than the Cholesky factor obtained from G(t) which is re- 
quired to be lower triangular (38|. Since the operators L{t) and p are diagonalizable 
in Fourier space, the matrix action of the operators U (t) and p on any vector can 
be computed using the Fast Fourier Transform with a cost of 0(N \og(N)). In sum- 
mary, our method allows in practice for the random variates of the stochastic driving 
field to be computed from g = Q(t)£ very efficiently, with a computational cost of 
only 0(N \og(N)). This is in contrast to the traditional Cholesky approach with a 
computational cost of 0(N 3 ). 

7. Regime II : Under-resolution of Hydrodynamic Relaxation (Quasi- 
Steady-State Limit). For many problems the equations of motion can be simplified 
by exploiting a separation of time-scales between the time-scale on which the hydrody- 
namic modes relax to a statistical steady-state and the time-scale associated with the 
motion of the microstructures. In this case the fluid equations can be approximated 

by 

(7.1) w=-L(i)- 1 [A + J]+a. 

The L = pLp T and the inverse is defined for the operator restricted to the linear 
space V = {w S R 3iV |pw = w}. The term a is introduced to account for the 
thermal fluctuations in this regime. We refer to this as the Quasi-Steady-State Stokes 
approximation [7J. Using this in equation 13. 3[ we obtain the following closed system 
of equations for the motion of the microstructures 

(7.2) =if SELM (t)[F]+J + A 
where 

(7.3) H SEhu (t) =-rZ(i)- 1 A 

(7.4) J = -rL(*)- x j. 

The A will be used to account for the thermal fluctuations. We consider the specific 
case when the operator A is linear in F and T is linear in u. In this case H s 

elm is a 

tensor which we refer to as the "effective hydrodynamic coupling tensor." 

In this regime the thermal fluctuations arise from the hydrodynamic modes which 
are relaxed to statistical steady-state. A key challenge is to determine the appropriate 
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statistics of A which accounts for the time integrated thermal fluctuations of the hy- 
drodynamics which impact the microstructure dynamics. For this purpose we rewrite 
equation 17.21 in differential form as 

(7.5) dX(i) = # SELM (t)Fdf + R(t)dB t 

neglecting for the moment J, and representing the contributions of A by R(t)dBf 
We derive the covariance structure S(t) — R{t)R{t) T by requiring consistency with 
the principle of Detailed-Balance of statistical mechanics [34[. The Fokker-Planck 
equation associated with equation 17.51 is 

(7.6) *jM~V.J 

(7.7) J = ff SELM (t)F* - ^S(t)Vx*. 

The ^(Xjt) is the probability density for the microstructures to have configuration 
X at time t. The equilibrium fluctuations of the system are required to have the 
Boltzmann distribution 

(7.8) *bd(X) = i cxp (-$(X)/fc fl T) 

where Z is a normalization constant which ensures the distribution integrates to 



34j |. Substituting this above and using that -Hselm is linear in F gives 



(7.9) J = (H SEhM (t) - 2j^r 5 (*)) F *bd 

where F = -Vx^- The principle of Detailed-Balance requires at thermodynamic 
equilibrium that J = 0. Requiring this to hold for all possible F gives 

(7.10) S(t) = 2k B TH SE ^ M {t). 

For S(t) to provide a covariance for a real-valued stochastic driving term, the hydro- 
dynamic coupling tensor H SEIjM (t) must be symmetric and positive semi- definite. In 
the case that A and T are linear operators this is ensured by condition [231 which from 
expression 17.31 gives 

(7.11) q T ff SELM (i)q = -v T (Z(i)- 1 ) vAx d > 0. 

To obtain this result we let v = T T q and use that L(i) -1 is symmetric negative 
definite. To obtain an approach useful in practice requires efficient methods for the 
generation of the stochastic driving term with covariance S(t). 

7.1. Generating the Stochastic Driving Field II. As discussed in Section|Hl 

a significant challenge in practice is to generate efficiently the Gaussian stochastic 
driving terms with the required covariance structure. We discuss an approach for 
SELM methods when the coupling operators A and T are linear. In this case 

(7.12) H SE ^ M (t) = -TL(t)- 1 T T Ax d 

by condition 13.41 Using properties of the operators discussed in Section [4J we can 
express the hydrodynamic coupling tensor as 

(7.13) ffs ELM (i) - (r(t)V(t)Ax d / 2 ) {r(t)V(t)Ax d / 2 Y . 
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We have used that the operators L(t) and p commute and since p is an exact pro- 
jection that p = p T , p = p 2 . Since L(t) is symmetric negative definite on the linear 
space V = {w e R 3Ar |pw = w}, we can factor L(t)" 1 = -V(t)V(t) T . The factor 
V(t) is readily obtained since L(t) and p are diagonalizable in the Fourier basis. From 
equations 17.131 and 17.101 we can factor the covariance as S(t) = R(t)R(t) T with 

(7.14) R(t) = (2k B TAx d ) 1/2 T(t)V(t). 

This expression for the factor can be used to compute the required Gaussian stochastic 
driving term A = R(i)£ with a computational cost of 0(N log(N) + M), where iV 
is the total number of mesh sites in the momentum field discretization and assuming 
the action of T can be computed with a cost of O(M) with M < N. 

The random variates are generated by utilizing the underlying discretization mesh 
of the momentum equations. This is accomplished by generating on the mesh uncor- 
rected standard Gaussian random variates £. Since V(t) is diagonal in the Fourier 
basis, the action V(t)£ is computed in Fourier space with a cost of only 0(N log(iV)). 
The operator T is then applied. If the operator T(t) makes use of only localized 
values of the mesh it can be computed with computational cost of 0{M). The last 
step in generating the random variate requires a scalar multiplication which incurs a 
computational cost of O(M). This procedure generates the stochastic driving term 
A with a computational cost of 0(Nlog(N) + M). For a sufficiently large number 
of microstructure degrees of freedom M, this method is significantly more efficient 
than the traditional approach based on Cholesky factorization of H 3EhM which costs 
0(M 3 ). 

7.1.1. Effective Hydro dynamic Coupling Tensor : H SELM . We now dis- 
cuss an approach for analyzing the effective hydrodynamic coupling tensors H SELM 
which appear in the quasi-steady-state formulation of the SELM approach. From 
equation 17.31 many types of hydrodynamic coupling tensors arc possible depending on 
the kinetic constraints and choice of coupling operators A and T. For concreteness we 
discuss the specific case corresponding to the Stochastic Immersed Boundary Method 
(SIB) ^HH- In the case of the SIB method, the specific coupling operators A and F 
are given by 13.101 and 13.111 From equation 17.31 the effective hydrodynamic coupling 
tensor is given by 

(7.15) 

[H IB (t) [F]]0=-X;a o (x m -x0(t)) 

m 

In the notation, the superscript [-]^ denotes for the composite vector the components 
associated with the j th microstructure degree of freedom. The [-] m denotes the vector 
components associated with the mesh site with index m. An analysis of variants of 
this tensor for point particles and slender bodies was carried-out in [J @| . 

Since H IB is linear in the microstructure forces, without loss of generality we 
can consider the case of only two microstructure degrees of freedom. We denote 
these as XM, X\ 2 ^ and the displacement vector by z = X\ 2 ^ — XW. In making 
comparisons with other hydrodynamic coupling tensors we find it helpful to make use 
of approximate symmetries satisfied by H IB . From equation 17. 151 H IB depends on 
z up to a shift of X^ relative to the nearest mesh site, and is similarly rotationally 



Lit)- 1 £)FMa a ( Xm -XM(t)) 
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Rotne-Prager-Yamakawa 
Effective IB Hydrodynamic Coupling Tensor Hydrodynamic Coupling Tensor 
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Fig. 7.1. Comparison of the Hydrodynamic Coupling Tensor of the Immersed Boundary Method 
H/b with the Oseen Tensor Hos and the Rotne-Prager-Yamakawa Tensor H n p Y . The components 
of the hydrodynamic coupling tensor for displacement r = |z| are shown for the parallel direction 
(circles) and the perpendicular direction (sguares). For two particles subject to an equal and opposite 
force, the velocity field corresponding to the Rotne-Prager- Yamakawa Tensor is shown on the right. 
The lighter shaded regions indicate a larger magnitude of the velocity. 



symmetry about the axis of z. This allows the tensor components for all configurations 
to be related to a canonical configuration with z = (z l7 0, 0). For any configuration 
this is accomplished by introducing the rotation matrix U so that Uz = (24,0,0) and 
considering H = UH IB U T . In our comparisons we consider H IB = (H), where the 
average is taken over all rotations and shifts with respect to the nearest mesh site. 

In practice, to numerically compute H IB we sample random configurations of 
XM and X™. A useful expression for the tensor components is iJy = ejHej = 
efv. In this notation, are the standard basis vectors in direction k and v is the 
microstructure velocity. For a computational implementation of the SELM method, 
this can be used by applying the force to the microstructure degrees of freedom 
and measuring the components of the realized microstructure velocities v. 

When using a SELM approach the hydrodynamic coupling tensor has features 
which depend on the discretization of the momentum equations, discretization of the 
microstructures, and the specific choice of coupling operators. For the specific choice 
of the IB coupling operators and discretization on a uniform mesh we discuss how the 
effective hydrodynamic coupling tensor compares with other hydrodynamic coupling 
tensors. We consider two specific tensors, the Oseen Tensor Q and the Rotne-Prager- 



Yamakawa Tensor [35J; |42|. The Oseen Tensor for a pair of particles experiencing 



equal and opposite forces can be expressed in terms of the displacement vector z as 

3 a 



H os (z) 



I ■ 



Similarly, the Rotne-Prager-Yamakawa Tensor can be expressed in terms of the dis- 
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placement vector z as 
2 

H RPY (z) = 

oirrja 

In the notation, 77 denotes the dynamic fluid viscosity, and a denotes the effective 
particle size in terms of the radius of a sphere. 

In Figure [7Tl the H IB is compared with the Oseen Tensor H os and Rotne-Prager- 
Yamakawa Tensor H RPY . It is found that the effective hydrodynamic coupling tensor 
of the Immersed Boundary Method agrees well with both of the tensors in the far-field 
r»d. An interesting finding is that in the near-field H IB shows very close agreement 
to H RPY , see inset in Figure ITU 

8. Applications. The SELM approach is expected to be applicable in the study 
of many different types of complex fluids and soft materials. As a demonstration of 
the proposed stochastic numerical methods, simulation studies are carried out for a 
few specific systems. These include studying: (i) the dependence of the shear viscosity 
on the shear rate in a FENE polymeric fluid, (ii) the frequency response of the elastic 
storage modulus and viscous loss modulus of a lipid vesicle fluid subject to oscillatory 
shear, (iii) the rheological responses over time of a gel-like material subject to a 
constant rate of shear. We now discuss each of these simulation studies in detail. 

8.1. Estimating Effective Macroscopic Stress. An important challenge in 
the study of complex fluids and soft materials is to relate bulk material properties 
to phenomena on the level of the microstructures of the material. To characterize 
properties of a material, experimental measurements are often made as a sample of 
material is subject to shear @; Q. To link microstructure mechanics, interactions, 
and kinetics to macroscopic material properties we develop estimators for an effective 
macroscopic stress tensor. Our estimators are based on similar approaches used to 
obtain the Irving-Kirkwood-Kramer formulas [5j;[6|;[l2; 25]. 

When using the SELM approach, the microstructures are modeled using n-body 
interactions and the domain is subject to generalized boundary conditions. For exam- 
ple, two body interactions can arise from bonds between monomer particles and three 
body interactions can arise from bond angle terms included in the potential energy 
Estimators for the stress must take these features into account. 

To obtain a notion of macroscopic stress we define a normal direction and a plane 
which cuts the unit cell. We then determine on average the forces exerted by the 
particles which lie above this plane on the particles which lie below this plane. We 
define the effective stress associated with this plane as the total of this exerted force 
divided by the area of the plane. To define an effective macroscopic stress we average 
over all possible planes within the unit cell having the specified normal direction, see 
Figure O 

More precisely, the effective macroscopic stress arising from n-body interactions 
is estimated using 




(8-1) <^L\J a A ^ (C)dC /' 

The L = b — a is the length of the domain in the z-direction and < ■ > denotes 
averaging over the ensemble. The A^™ denotes the microscopic stress arising from 
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Fig. 8.1. An effective macroscopic stress is estimated from a sample of the material by com- 
puting the forces transmitted across a plane which cuts through the sample at a specified location 
and with a specified normal. At the level of the micro structures, the cut-plane is used to divide the 
sample into two bodies labeled A and B, shown in the middle. The effective stress is estimated by 
computing the force exerted by particles in body A on particles in body B. For models with two body 
interactions a contribution is made to the stress only if one particle is in body A while the other is 
in body B, shown on the far right on top. For three body interactions there are two possible cases 
for how forces can be transmitted across the cut plane, shown on the far right in the middle and 
bottom. 



the n-body interactions associated with a given stress plane and is defined by 

n— 1 k k 



1 n— 1 k k n 



qgQ„fc=lj=l j=l j=k+l 

The Q n is the set of n-tuple indices q = (qi, . . . ,q n ) describing the n-body interactions 
of the system, f q j denotes the force acting on the j th particle of the interaction, and 
x 9j denotes the j th particle involved in the interaction. As a matter of convention in 

(z) (z) 

the indexing q we require that i < j implies x^ < x. qj . This expression corresponds 
to a sum over all the forces exerted by particles of the material above the cross-section 
at £ — z on the particles of the material below. Each term of the summation over 
k = 1, . . . , n— 1 corresponds to a specific number of particles of the n-body interaction 
lying below the cross-section at £ = z, see Figure [QI 

When integrating the microscopic stress, a useful identity is that 



f n* =1 w(c - xW) • ny =fc+1 w(xW - CR = x* ; « - x*f > 

J a 



(8.3) 
where 



(8.4) x*f ) = 

By integrating equation 18.21 we obtain 

(^) fASlm^EEEfi't^ 

Ja q eQ„fc=ij=i 
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Fig. 8.2. Components of the second moment of the extension vector are shown as the shear 
rate is varied. The second moment matrix is M = (zz T ). On the left is shown the off diagonal 
entry Mi,3 as a function of shear rate. On the right is shown the averaged mean squared extension 
vector of the dimer, which is given by P = (|z| 2 ) = Trace[M]. The moments show a significant 
dependence on the rate of shear. 



This can be further simplified by switching the order of summation of j and k and 
using the telescoping property of the summation over k. From equation l8 . 21 this yields 
the following estimator for the stress contributions of the n-body interactions 

(8-6) ^4EE(«'(^-f))- 

qeg„ j=i 

This defines an effective macroscopic stress tensor contribution in terms of the n-body 
interactions of the microstructures of the material. To obtain the total contribution 
of the microstructure interactions to the stress, all of the contributions of the n-body 
interactions are summed to obtain the effective macroscopic stress tensor 

(8.7) ^=e^- 

n 

This notion of the macroscopic stress will be used to link bulk rheological properties 
to the microscopic simulations. 

8.2. Application I: Complex Fluid of Finite Extensible Non-linear Elas- 
tic (FENE) Dimers. As a demonstration of the proposed computational method- 
ology we consider a fluid with microstructures consisting of elastic polymers. The 
polymers are modeled as idealized elastic dimers which have the potential energy 

(8.8) 4>{r) = \kt% log^l- (j^ ^ . 

The K denotes the polymer stiffness, r denotes the length of extension of the dimer, 
and ro denotes the maximum permitted extension length The configuration of 
each dimer will be represented using two degrees of freedom , X( 2 ) . The potential 
energy for the dimer is given by ^(X) = 0(|X( 2 ) — X^l), where X is the composite 
vector for the particle configuration. 
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Parameter 


Description 


N 


Number of mesh points in each direction. 


Ax 


Mesh spacing. 


L 


Domain size in each direction. 


T 


Temperature. 


k B 


Boltzmann's constant. 


P 


Dynamic viscosity of the solvent fluid. 


P 


Mass density of the solvent fluid. 


K 


Bond stiffness. 


TO 


Maximum permissible bond extension. 


Is 


Stokesian drag of a particle. 


7° 


Shear rate amplitude. 


7° 


Strain rate amplitude. 


a 


Effective radius of particle estimated via Stokes drag. 



Table 8.1 

Description of the parameters used in simulations of the FENE polymeric fluid. 



Parameter 


Value 


TV 


36 






Ax 


11.25 nm 






L 


405 nm 






T 


300 K 






k B 


8.3145 x 10 3 


nm 2 


amu • ns -2 • K" 1 


P 


6.0221 x 10 5 


amu 


cm -1 • ns -1 


P 


6.0221 x 10 2 


amu 


nm 


K 


8.9796 x 10 3 


amu 


ns^ 2 


ro 


200 nm 






Is 


1.7027 x 10 8 


amu 


ns" 1 


a 


15 nm 







Table 8.2 



Values of the parameters used in simulations of the FENE polymeric fluid. 



When the polymeric fluid is subject to shear the thermally fluctuating polymeric 
microstructures are expected to significantly re-orient and deform as a consequence 
of the shear stresses. This along with thermal fluctuations of the microstructures is 
expected to play an important role in the bulk response of the polymeric fluid. To link 
the bulk material properties of the fluid to the microstructures, we use the effective 
macroscopic stress cr p obtained from equation l8.6l To characterize the bulk rheological 
response we consider the shear viscosity rj„ and first normal stress coefficient of 
the polymeric fluid. We define these as 0] 

(8-9) fb = 4'' V) li 

(8.10) = tr^/T 2 . 

The 7 is the rate of shear of the polymeric fluid. In the notation, the superscript 
(s,v) indicates the tensor component with the index s corresponding to the direction 
of shear and the index v corresponding to the direction of the fluid velocity. The 
contributions of the solvent fluid to the shear viscosity and normal stresses can be 
considered separately Q. 
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Shear Viscosity First Normal Stress Coefficient 




Fig. 8.3. Rheological properties of the FENE polymeric fluid are shown as the rate of shear is 
varied. The shear viscosity is shown on the left and the first normal stress difference is shown on 
the right. As the shear rate increases the dimers align increasingly with the direction of fluid flow, 
shown as insets. 



The SELM approach is used to study how the shear viscosity and first normal 
stress difference depend on the rate of shear of the polymeric fluid. Simulations are 
performed using the SELM method in the regime where the hydrodynamic modes 
are relaxed to statistical steady-state with parameters given in Table 18.11 For A 
and r the coupling tensors of equation 13.101 and equation 13.111 are used. From an 
ensemble average over many computational experiments the moments of the extension 
vector z are estimated as the shear rate is increased. The polymeric microstructure 
moments are seen to respond strongly as the shear stresses of the fluid increase, see 
Figure 18.21 This indicates that the rheological properties of the polymeric fluid will 
depend significantly on the rate of shear. The SELM simulations show that the shear 
viscosity and the first normal stress difference do in fact vary significantly with the 
shear rate, see Figure IQ1 

The shear viscosity is found to decrease as the shear rate increases. This appears 
to occur as a consequence of the dimers increasingly aligning with the direction of 
the fluid flow and as a consequence of the dimers approaching the maximal extension 
permitted by eauation l8.8l The increased extension results in a non- linear increase in 
the effective stiffness of the dimer (defined for a given extension by Taylor expanding 
to second order equation 18. 8p . While the dimers become increasingly extended with 
stronger restoring forces this is counter-balanced by the dimers being increasingly 
stiff and the thermal fluctuations less frequently driving the dimer into configurations 
crossing the stress plane. The net effect is that the mechanical stress transmitted on 
average by the dimers in the direction of shear does not increase as the shear rate 
increases. This results in a lower effective shear viscosity (note the division by 7 in 
equation 18. 9p . This is a well-known phenomena in polymeric fluids and is referred to 
as shear thinning. The simulations demonstrate that the SELM approach is capable 
of capturing at the level of the microstructures such phenomena, see Figure 18.31 

8.3. Application II: Polymerized Lipid Vesicle Fluid. As a further demon- 
stration of the applicability of the SELM approach we show how the stochastic nu- 
merical methods can be used to investigate the bulk material properties of a complex 
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Fig. 8.4. Recursive Method for Mesh Construction. The triangulated mesh for a spherical 
vesicle is constructed by starting with the vertices and faces of a regular icosahedron, shown on the 
left. The edges of the icosahedron are bisected and connected to divide each triangular face into four 
smaller triangular faces. The vertices located at the bisection points are projected radially outward 
to the surface of the sphere, shown in the middle. This refinement procedure is repeated recursively 
until a mesh of sufficient resolution is obtained. The mesh obtained after two levels of recursive 
refinement, which we use to represent polymerized vesicles, is shown on the right. 

fluid with polymerized vesicle micro-structures. We discuss how the methods can be 
used to compute the response of the complex fluid subject to an oscillating shear flow 
varied over a wide range of frequencies. 

To obtain a triangulated mesh which captures the shape of a vesicle having a 
spherical geometry we start with an icosahedral which is circumscribed by a sphere 
of a given radius. We use the faces of the icosahedron as an initial triangulated mesh. 
To obtain a mesh which better approximates the sphere we bisect the three edges of 
each triangular face to obtain four sub-triangles. The newly introduced vertices are 
projected radially outward to the surface of the sphere. The process is then repeated 
recursively to obtain further refinements of the mesh. This yields a high quality mesh 
for spherical geometries. A vesicle represented by a mesh obtained using two levels of 
recursive refinement is shown in Figure WM 

To account for the mechanics of a polymerized vesicle the following interactions 
are used for the control points of the mesh 



The r denotes the displacement between two control points, I denotes a preferred 
distance between control points, and r denotes a normalized displacement vector 
(tangent vector) between two control points. The <f>\ energy accounts for the stretching 
of a bond between two control points beyond its preferred extension. The 4>2 energy 
accounts for bending of the surface locally by penalizing the misalignment of tangent 
vectors. 

For a given triangulated mesh of control points the total energy is given by 



(8.11) 



(8.12) 



<Mti,T 2 ) = o^l l T l - T : 




(8.13) 
(8.14) 





(<,3')eQi 



(8.15) 





(i,j,fc)eQ 2 



22 



P.J. ATZBERGER 



Low Frequency Response 




= 1.6* = O.O?r 6 = 0.4jt 



2 



1 













a 



•1 



■2 







0.5 



1 

eft 



1.5 



2 



Fig. 8.5. Simulation results showing the vesicle response when subject to an oscillating shear 
flow. At low frequency the vesicle shape distortion is small and is masked by thermal fluctuations. At 
low frequency the vesicle membrane stresses equilibrate to a good approximation with the bulk shear 
stresses, as illustrated in the plot of a xz . For the vesicle configurations shown, the low frequency 
response corresponds to u> = 3.9294 X 10~ 3 ns~ 1 , 7 = 1.9647 X 10 -3 ns -1 , <ro = 3.7114 X 10 s amu- 
nm~ x ■ ns~ 2 . The phase 9 = u>t is reported in the range [0, 2n). For additional parameters used in 
the simulations see Table HOI anri fg^l 

The X denotes the composite vector of control points. The j th control point is denoted 
by X^l . The Qi and Q2 are index sets defined by the topology of the triangulated 
mesh. 

The first energy term E\ accounts for stretching of the vesicle surface and is 
computed by summing over all local two body interactions Q± defined by the topology 
of the triangulated mesh. For the distance = |X^ — X\^ \ between the two points 
having index i and j, the energy E\ penalizes deviations from the preferred distance 
£ij. The preferred distances tij are defined by the geometry of a spherical reference 
configuration for the vesicle. To ensure the two body interactions are represented by 
a unique index in Qi we adopt the convention that i < j. 

The second energy term E2 accounts for curvature of the vesicle surface and 
is computed by summing over all local three body interactions Q2 defined by the 
topology of the triangulated mesh. The energy penalizes the the misalignment of the 
tangent vectors = (X^ - XM)/^ and T jk = (X^ - X^)/r jk . In the set of 
indices in Q2 it is assumed that the point with index j is always adjacent to both i 
and k. To ensure the three body interactions are represented by a unique index in Q2 
we adopt the convention that i < k. 

To investigate the bulk rheological properties, the complex vesicle fluid is sub- 
jected to an oscillatory shear with rate 7 = 7° cos(wt). We consider the dilute regime 
in which it is sufficient to study a single polymerized vesicle subject to oscillatory 
shear. To estimate the effective macroscopic stress tensor the tensor is decomposed 
into contributions from two body and three body interactions 



For the contributions of the n-body interactions to the macroscopic stress a\ J we use 



(8.16) 
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High Frequency Response 




Fig. 8.6. Simulations results showing the vesicle response when subject to an oscillating shear 
flow. At high frequency the vesicle shape is visibly distorted and the membrane stresses do not 
have time to equilibrate with the bulk shear stresses, as illustrated by the configurations for phase 
9 = 1.6,0.4 and the plot of cr xz . For the vesicle configurations shown, the high frequency response 
corresponds to u) = 1.2426 X 10 2 ns _1 , 7 = 6.2129 X 10 1 ns -1 , o = 4.6314 X 10 10 amu - nm~ x ■ ns~ 2 . 
The phase = u)t is reported in the range [0, 2tt). For additional parameters used in the simulations 
see Table \8j\and fO 



the approach discussed in Section 18.11 and the specific estimator given by equation 

m 

For many materials, the responses of the stress component a xz (t) to bulk stresses 
and strains are linear to a good approximation over a wide range of frequencies 



provided the stress and strain amplitudes are sufficiently small 32]. As a mea- 
sure of the material response, we consider the dynamic complex modulus G(ui) = 
G'(oj) + iG"(u>), whose components are defined from measurements of the stress as 
the best least-squares fit of the periodic stress component cr xz (t) by the function 
g(t) = G'(oj)j° cos(wt) + G"(uj)j° sin(wt). This offers one characterization of the re- 
sponse of the material to oscillating bulk shear stresses and strains as the frequency 
ui is varied. 

To estimate the dynamic complex modulus in practice the least squares fit is 
performed for a xz (t) over the entire stochastic trajectory of the simulations (after 
some transient period). Throughout our discussion we refer to 6 = Lut as the phase 
of the periodic response. In our simulations the maximum strain each period was 
chosen to always be half the periodic unit cell in the x-direction, corresponding to 
strain amplitude 7 — i. This was achieved by adjusting the shear rate amplitude 
for each frequency using 7 = 7°cj. 

Simulations were performed with the SELM approach in the regime where the 
hydrodynamic modes were treated as relaxed to statistical steady-state. The specific 
coupling operators A and T from l3.10l and l3.lT1 were used. The simulation results of the 
complex modulus response of the vesicle when subject to a wide range of frequencies is 
shown in Figure EH Figurc lSlfl and Figure E3 It was found that at low frequency the 
vesicle shape distortion is small and masked by thermal fluctuations. At low frequency 
the vesicle membrane stresses equilibrate to a good approximation with the bulk shear 
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Elastic Storage Modulus Viscous Loss Modulus 




10"" 10" 2 10° 10 z 10* 10"* 10" 2 10° 10 2 10* 



Fig. 8.7. Frequency response of the dynamic complex modulus of the vesicle fluid subject to 

an oscillating shear flow. Throughout the simulations the total strain was held fixed to be half the 

domain length, 7 = \L. For a description of the parameters and values used in the simulations, 
see TableUnH and[8^\ 



stresses, as illustrated in the plot of a xz in Figure [53] It was found at high frequency 
the vesicle shape is visibly distorted and the membrane stresses do not have time to 
equilibrate with the bulk shear stresses, as illustrated by the configurations for phase 
9 = 1.6, 0.4 and the plot of a xz in Figure [831 For the vesicle configurations shown, the 
low frequency response corresponds to u> — 3.9294x 10 _3 ns _1 , 7 = 1.9647x 10 _3 ns _1 , 
(To = 3.7114 x 10 8 amu • nm^ 1 • ns~ 2 , and the high frequency response corresponds to 
lu = 1.2426 x 10 2 ns-\ 7 = 6.2129 x lC^ns^ 1 , a = 4.6314 x 10 10 amu ■ nm" 1 ■ ns" 2 . 
The phase 9 = uut is reported in the range [0, 2tt). A description of the parameters 
and specific values used in the simulations can be found in Table 18.31 and 18.41 



Parameter 


Description 


N 


Number of mesh points in each direction. 


Ax 


Mesh spacing. 


L 


Domain size in each direction. 


T 


Temperature. 


k B 


Boltzmann's constant. 




Dynamic viscosity of the solvent fluid. 


P 


Mass density of the solvent fluid. 


K x 


Vesicle bond stiffness. 


K 2 


Vesicle bending stiffness. 


D 


Vesicle diameter. 


LU 


Frequency of oscillating shearing motion. 


9 


Phase of the oscillatory motion, 9 = tut. 


7 


Shear rate. 


7° 


Shear rate amplitude. 


7 


Strain rate. 


7° 


Strain rate amplitude. 



Table 8.3 

Description of the parameters used in simulations of the vesicle fluid. 
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Parameter 


Value 


N 


27 


Ax 


7.5 nm 


L 


2.025 x 10 2 nm 


T 


300 K 


k B 


8.3145 x 10 3 nm 2 • amu • ns~ 2 • KT 1 


i 1 


6.0221 x 10 5 amu • cm" 1 • ns" 1 


P 


6.0221 x 10 2 amu • nm" 3 


Ki 


2.2449 x 10 7 amu • ns~ 2 


K 2 


8.9796 x 10 7 


D 


50 nm 



Table 8.4 



Fixed values of the parameters used in simulations of the vesicle fluid. 
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Fig. 8.8. Model for a Gel-Like Material. The gel is formed by polymeric chains which bind 
together, shown on the left. The polymeric chains are each comprised of five control points and 
have specialized binding sites at the second and fourth control point, shown in the center. The inter- 
polymer bonds have a preferred extension and angle. When an inter-polymer bond is strained beyond 
50% of its preferred rest length the bond breaks irreversibly, shown on the right. 



8.4. Application III: Rheology of a Gel-like Material. As a further demon- 
stration of the applicability of the SELM approach we show how the stochastic nu- 
merical methods can be used to investigate properties of a gel-like material subject 
to shear. The methods are used to study how the shear viscosity changes over time 
as the gel is subjected to shear at a constant rate. 

The gel-like material is modeled as a collection of polymer chains which are able 
to bond together at two specialized sites along the chain, see Figure 18.81 The energy 
associated with the mechanics of the individual polymer chains and the bonds which 
they form are given by 



.17) 
.18) 

.19) 
.20) 



h(ri,T 2 ) = -K 2 |ti - r 2 | 2 



2 

h(r) = a 2 K 3 exp 
b 4 {9) = -K 4 cos(t 



(r - r 0i3 ) 2 



2a 2 



7 0,4J 



The r is the separation distance between two control points, 9 is the bond angle 
between three control points, and r is a tangent vector along the polymer chain, see 
Figure [87H1 

The 4>\ energy accounts for stretching of a bond within a polymer chain from its 
preferred extension ro.i. The </>2 energy accounts for bending of the polymer chain 
locally. To account for interactions at the specialized binding sites of the polymers 
the potentials 03 and 04 are introduced. The potential 03 gives the energy of the 
bond between the two polymer chains and penalizes deviation from the preferred 
bond extension tq^. The exponential of 03 is introduced so that the resistance in the 
bond behaves initially like a harmonic bond but decays rapidly to zero when the bond 
is stretched beyond the length a. The potential 04 gives the energy for the preferred 
bond angle when two of the polymer chains are bound together. 

The total energy of the system is given by 

(8.21) $[X] = £ : [X] + E 2 [X] + E 3 [X] + E 4 [X\ 
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Shear Viscosity vs Time 




5000 10000 
time (ns) 



Fig. 8.9. Study of the shear viscosity of a gel-like material. At time zero the material has weak 
bonds between short polymeric chains. Under the shear deformation the gel is stretched and the bonds 
are strained until ultimately breaking. Many of the polymers are misaligned with the direction of 
fluid flow and are further stretched by the fluid shear stresses. As the polymer chains align with the 
direction of fluid flow the forces transmitted in the direction of shear decrease and the shear viscosity 
approaches a steady-state value. The thermal fluctuations maintain transient misalignments of the 
polymer chains which transmit forces in the direction of shear resulting in a contribution to the 
shear viscosity which is non-zero at steady-state. The micro structure reordering in each of these 
stages, labeled I, II, III, is reflected in the shear viscosity of the material as a function of time and 
in Figure Y8.10\ For the specific physical parameters used in these simulation see Table l8^5\ and l876\ 



.22) e 1 [x}= J2 Mnj), e 2 [x} = Yl Hra.Tjk) 

(ij)eSi (ij,fe)eQ 2 
.23) E 3 [X} = J2 &(r«), E 4 [X}= £ 



The sets Qk define the interactions according to the structure of the individual poly- 
mer chains and the topology of the gel network. When bonds are stretched beyond 
the critical length 3<t they are broken irreversibly, which results in the sets Q3 and 
Qi being time dependent. 

To study the rheological response of the gel-like material the system is subjected 
to shear at a constant rate. To obtain an effective macroscopic stress a v for the system 
the estimator is used from equation I8.6I To characterize the rheological response we 
use the shear viscosity defined by 



(8.24) 



ft- 



The 7 is the rate of shear of the polymeric fluid. In the notation, the superscript (s, v) 
indicates the tensor component with the index s corresponding to the direction of shear 
and the index v corresponding to the direction of the fluid velocity. The contributions 
of the solvent fluid to the shear viscosity can be considered separately Q . 

The entire gel network experiences an unbounded shear deformation. This is 
expected to result in breakage of bonds of the gel network. This suggests that the 
rheological response will depend on how long the material has been subject to shear. 
To investigate the role reorganization at the microstructure level, repeated stochastic 
simulations are carried out using the SELM approach to determine the effective shear 
viscosity of the material as a function of time. 
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Fig. 8.10. The micro structure of a gel-like material at three different times. On the left is shown 
the micro structure of the gel-like material before any shear has been applied. In the middle is show 
the micro structure of the gel after almost all of the bonds between polymer chains have been broken. 
In this case, the misaligned polymer chains continue to be stretched by the shear stresses of the fluid 
yielding a relatively large effective shear viscosity. On the right is shown the micro structure of the 
gel when the system has relaxed to statistical steady-state. In this case, the thermal fluctuations 
drive transient misalignments of the polymer chains with the direction of flow which on average 
make a non-zero contribution to the shear viscosity. The times shown in each of these figures is 
t = ns, t = 2844 ns, t = 7111 ns. For the specific physical parameters used in these simulations 
see Table[8li\ and[8~6\ 



Parameter 


Description 


N 


Number of mesh points in each direction. 


Ax 


Mesh spacing. 


At 


Time step. 


L 


Domain size in each direction. 


T 


Temperature. 


k B 


Boltzmann's constant. 




Dynamic viscosity of the solvent fluid. 


P 


Mass density of the solvent fluid. 


7 


Shear rate. 


N p 


Number of polymer chains. 


N s 


Number of control points per polymer chain. 


r p 


Polymer effective cylindrical radius. 


Ki 


Stiffness of the bonds of the polymer chain. 


r 0,l 


Rest length of the bonds of the polymer chain. 


K 2 


Bending stiffness of the polymer chain. 


K 3 


Stiffness of the bonds at a polymer binding site. 


m,3 


Rest length of the bond at a polymer binding site. 


K 4 


Bending stiffness of the bond at a polymer binding site. 


#0,4 


Preferred angle of a bond at a polymer binding site. 



Table 8.5 

Description of the parameters used in simulations of the gel-like material. 



An interesting behavior is found in which the material initially exhibits an in- 
creased shear viscosity before settling down to a steady-state value. The responses of 
the material to shear can be roughly divided into three stages. In the first, there is 
an initial increase which can be attributed to the stretching of the inter-chain bonds 
between the polymer chains and the intra-chain bonds within each polymer chain, 
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Parameter 


Value 


TV 


72 






Aa; 


11.25 nm 




At 


1.4222 ns 




L 


810 nm 






T 


300 K 






k B 


8.3145 x 


10 3 nm 2 


amu • ns -2 • K -1 


A* 


6.0221 x 


10 5 amu 


cm -1 • ns -1 


P 


6.0221 x 


10 2 amu 


nm -3 


7 


1.2 x 10" 


3 ns" 1 




N p 


110 






N s 


5 






r P 


15 nm 






K X 


2.9932 x 


10 5 amu 


ns -2 


ro,i 


30 nm 






K 2 


2.9932 x 


10 8 




K 3 


2.9932 x 


10 5 amu 


ns z 


ro.3 


30 nm 






K 4 


2.9932 x 


10 8 




#0,4 


70° 







Table 8.6 



Fixed values of the parameters used in simulations of the gel-like material. 



which occurs as the gel as a whole is strained. After a relatively short period, the 
bonds between the polymer chains are observed to break with the remaining stress 
arising from the stretching of the polymer chains which occurs from the shear stresses 
of the fluid and misalignment with the direction of flow, see the region labeled by I 
in Figure \SM and EM 

In the second stage, the individual polymer chains rotate and begin to align with 
the direction of flow. As a result of the intra-chain restoring forces the strain of the 
individual polymer chains is reduced. The increased alignment and reduced strain of 
the polymer chains yields an overall decrease in the forces transmitted in the direction 
of shear. Consequently, the shear viscosity begins to decrease, see the region labeled 
by II in Figure [EU and EHUl 

In the last stage, the chains eventually settle into a statistical steady-state in 
which the thermal fluctuations drive the chains to misalign only transiently with the 
flow direction. These misaligned excursions by the polymer chains sustained by the 
thermal fluctuations result in forces transmitted in the direction of shear on average. 
This is reflected in the shear viscosity by a non-zero steady-state value, see the region 
labeled by III in Figure [EH and EM 

Using the SELM approach more complicated situations could also be studied, 
such as the case in which the bonds between the polymer chains are able to reform. 
An interesting investigation in this case would be to study how the viscosity behaves 
after decreasing or ceasing shearing of the system for a period of time. In this case the 
gel would have time to reform structures before being again subjected to large shears. 
Using such a SELM approach a widely variety of shear thinning and thixotropic 
phenomena could be studied at the level of the microstructures [1; H; [Hj] . 
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9. Conclusions. A general formalism was developed which allows for the cou- 
pling of Eulcrian and Lagrangian descriptions of physical systems. A general approach 
was introduced for incorporating thermal fluctuations in such descriptions. The ap- 
proach addresses both the inertial regime and the overdamped regime. For the study 
of rheological responses of materials, an approach was developed which allows for gen- 
eralized periodic boundary conditions which induce the shear. For simulations using 
the formalism stochastic numerical methods were developed which efficiently generate 
the required stochastic driving fields. As a demonstration of how these methods can 
be used in practice, simulation studies were carried out for complex fluids and soft 
materials. The basic Stochastic Eulerian Lagrangian Method (SELM) approach is ex- 
pected to be useful in the formulation of descriptions and computational approaches 
for the study of a wide variety of fluid structure phenomena involving thermal fluctu- 
ations. 
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Appendix A. Invariance of the Boltzmann Distribution under SELM 
Stochastic Dynamics. 

The probability distribution of the stochastic equations 12.1112.21 are governed for- 
mally by the Fokker-Planck equation 

(A.1) - = - V.J 



with the probability flux given by 
(A.2) J = 



£p* + (A + A)*-±GV p * 
(r + 7 )#- iVKVx* 



The G, W are the covariance operators associated with g and Z. The fy(p,~K,t) is 
the formal probability density for finding the system in state (p, X) at time t. For 
the present purposes our discussion will only be formal since the SPDEs are infinite 
dimensional and for the density there is no Lebesgue measure for the function space, 
see 11; ]1| 3(J. In practice a finite dimensional stochastic process will always be used 
to approximate the SPDEs and has a probability distribution satisfying a well-defined 
equation. 

For the systems under consideration, the Boltzmann distribution has the form 
^ B£>(P)X) = i exp [— E\p, X]/fe_BT], where Z is a normalization constant so that 
^ bd integrates to one [34| . The requirement that this distribution is invariant under 
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the stochastic dynamics of !2.Hf^TS1 is equivalent to V ■ J = 0. This requires 



(A.3) 
(A.4) 
(A.5) 



V • J 



Ax 
A 2 



Ax + A 2 + V • A 3 = 

[(A + A) • V P E + (r + 7 ) • V x £] (-fcsT)- 1 * 
(V p -(A + A) + Vx-(r + 7 ))* 



(A.6) 



A, 



■3 




For the energy given 



by equation 12.31 we have 



(A.7) 
(A.8) 



V p E = Po x p 
V X E = V x $ 



F 



where F denotes the force for the configuration. 

Now we can derive conditions for the coupling operators by requiring that Ax = 
A 2 = for all possible values of p and F. The requirement that Ax = corresponds to 
the energy being conserved under the dynamics of equations 12.1112.21 when g = Z = 
and er = 0. For these dynamics the energy satisfies dE/dt = (A + A) • W p E + (T + 
7) • Vx-E = 0. Since the forces associated with time independent constraints do not 
do any work on the system we have that A • V P E + 7 • Vx-E = 0. Conservation of 
energy then requires A • V p i? + T ■ Vx-E = 0. By using the variational derivatives 
of E given in fOTOl we have A ■ V p E = J Ap^pdx and T ■ V x £ = / -TFdq. By 
substituting these expressions into IA. 41 we obtain from Ax = that the condition l2.4l 
must be satisfied. 

The requirement that A 2 — requires that the dynamical flow in phase space 
defined by (A + A, T + 7) is volume preserving. For the dynamics when g = Z = 
0, <t = 0, and E = this condition is equivalent to requiring that the uniform 
distribution is invariant under the dynamics. The condition 12.61 follows by using 
the function representing the variational derivatives (20| appearing in the divergence 
operation corresponds to Vx ■ T — J*(£T/5X)(q, q)dq, Vx ■ 7 = /(<57/<5X)(q, q)c£q, 
and similarly for A, A. 

The requirement that A3 — requires from equation IA.7IIA.8l that 
£p + [(Gp~ 1 p — WF)/2kBT] — for any p and F. This requirement corresponds 
to the condition of Detailed-Balance of statistical mechanics [34| . Since p and F are 
arbitrary, this requires that W — so that Z = 0. This also requires that G = —2CC 
with C = ksTpoT, where X is the identity operator. This yields condition 12.71 From 
the form of the energy in 12.31 and the Boltzmann distribution we see the equilibrium 
fluctuations of p are Gaussian with covariance C. This condition relates the equi- 
librium fluctuations to the dissipative operator of the system and is a variant of the 
Fluctuation-Dissipation Principle of statistical mechanics 34]. This shows that pro- 
vided the coupling operators and stochastic fields satisfy conditions 12.41 12.61 and 12.71 
the Boltzmann distribution is invariant under the SELM stochastic dynamics. 

For the discretized equations, we now derive conditions 13.41 13.61 and 13.81 The 
calculations follow similarly to the case above so we only state the basic features of 
the derivation. For the discretized equations the probability flux is given by 



(A.9) 



J 



Lp* + (A + A)* 

(r + 7 )* 
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where p and X are now finite dimensional vectors. The Boltzmann distribution now 
uses the energy of the discrete system 

(A.IO) E[p, X] = £ lpo 1 \Pm\ 2 Ax d + $(X) 

m 

with 

(A.ll) X7 p E = p 1 pAx d 

(A.12) V x £ = Vx$ = -F. 

Substituting these expressions in IA.3I - IA.6I and reasoning as above yields the condi- 
tions E31 EH andEU 

Appendix B. A Fluctuation-Dissipation Principle for Time-Dependent 
Operators. Consider the stochastic process given by 



(B.l) dz t = L(t)zdt + Q{t)dB t 

(B.2) G(t) = QQ T ■ 

We now establish the following fluctuation-dissipation relation 

(B.3) G(t) = -L(t)C - C T L(t) T . 



This relates the covariance G(t) of the stochastic driving field to a time-dependent 
dissipative operator L{t) and a time-independent equilibrium covariance C . We show 
that this relation allows for G(t) to be chosen to ensure that the stochastic dynamics 
exhibits at statistical steady-state equilibrium fluctuations with the specified covari- 
ance C. 

Let the covariance at time t be denoted by 

(B.4) C(t) = (u(t)u(tf). 

By Ito's Lemma the second moment satisfies 

(B.5) dC(t) = (L(t)C(t) + C(t) T L(t) T + G(t)) dt. 

It will be convenient to express this equation by considering all of the individual entries 
of the matrix C(t) collected into a single column vector denoted by c t . Similarly, for 
covariance matrix G(t) we denote the column vector of entries by gt and for C by c. 
Since the products L(t)C(t) and C(t) T L(t) T are both linear operations in the entries 
of the matrix C(t) we can express this in terms of multiplication by of a matrix A(t) 
acting on c t . 

This notation allows for equation IB. 51 to be expressed equivalently as 
(B.6) dc* = (A(t)c t + g t ) dt. 

The equation lB.5l can be solved formally by the method of integrating factors to obtain 

(B.7) ct =e H (°^c + [ e s ^ Ss ds 

Jo 

where H(s, t) — f A(r)dr. 
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The fluctuation-dissipation relation given by equation IB. 31 is equivalent to choos- 
ing 

(B.8) gs - -A(s)c. 

For this choice, a useful identity is 

(B.9) e s ^g s = ^-e s ^c. 

OS 

Substitution into eauation lB.7l gives 

(B.IO) c t = e s (°'')c + (e 5 ^ - e H ( ^) c. 

Now, if L(t) is negative definite uniformly in time, v T L(t)v < ao < 0, then A(t) 
is also uniformly negative definite. This implies that 

(B.ll) lim e H (°'*) = 0. 

t — >oo 

Taking the limit of both sides of equation IB. 101 and using equation IB. Ill yields 
(B.12) limc t =c. 

t — *oo 



This shows that the stochastic driving field with covariance given by equation IB. 31 
yields equilibrium fluctuations with covariance C . This extends the fluctuation- 
dissipation relation to the case of time-dependent operators. 

For the discretization given in Section [4j we point out some of the properties of 
the specific matrix L(t) which are used. From equation 14. 121 the non-zero eigenvalues 
of L(t) can be shown to be negative and uniformly bounded away from zero in time. 
The eigenvector associated with the zero eigenvalue of L(t) is in fact the same for 
all times. The eigenvector of the zero eigenvalue is proportional to the vector with 
all components set to one. In practice, this mode is set to zero. By conservation 
of momentum of the fluid body as a whole, this mode remains zero when subject to 
internal conservative forces. This allows for the operator L(t) to be considered as 
acting on the linear space which excludes this null eigenvector. On this linear space, 
L(t) is strictly negative definite uniformly in time. Similar considerations can be made 
when considering the effect of the incompressibility constraint for the operator L(t) = 
pL{t). Thus the time-dependent fluctuation-dissipation relation given by equation lB.3l 
still holds provided the stochastic process is considered on the appropriate linear space 
which excludes the null eigenvectors. 

Appendix C. The Particle Representation Function S a . In the immersed 
boundary method, it is required that a function 5 a be specified to represent the 
elementary particles. The representation of this function is often derived from the 
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following function (f> which is known to have desirable numerical properties 0; l3~fl | : 



r 



(C.l) 



if r < -2 



| (5 + 2r - V-7 - 12r - 4r 2 ) , if -2 < r < -1 

, if -1< r < 



| (3 + 2r + Vl -4r-4r 2 ) 
§ (3 - 2r + Vl + 4r - 4r 2 ) 



if < r < 1 



| (5 - 2r - V-7+12r-4r 2 ) , if 1 < r < 2 



I 



if 2 < r. 



For three dimensional systems the function 5 a representing elementary particles 
of size a is 



(C.2) 



1 



.(2) 



.(3) 



where the superscript indicates the index of the vector component. 

To maintain good numerical properties, the particles are restricted to sizes a = 
nAx, where n is a positive integer. For a derivation and a detailed discussion of the 
properties of these functions see [1; 31 1 . 



Appendix D. Table. 



Parameter 


Description 


N A 


Avogadro's number. 


amu 


Atomic mass unit. 


nm 


Nanometer. 


ns 


Nanosecond. 


Ub 


Boltzmann's Constant. 


T 


Temperature. 


V 


Dynamic viscosity of water. 




Stokes' drag of a spherical particle. 






Parameter 


Value 




N A 


6.02214199 x 1(F. 




amu 


1/10 3 AU kg. 




nm 


1(T 9 m. 




ns 


1CT 9 s. 




k B 


8.31447 x 10 3 amu nm 2 /ns 2 K. 




T 


300K. 




T) 


6.02214199 amu/cm ns. 





